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ABSTRACT 


Numerical  studies  of  time-dependent,  two-dimensional,  stratified 
flow  of  incompressible,  viscous,  diffusive  fluid  are  described.  The 
particular  physical  problem  which  has  been  considered  is  that  of  strati¬ 
fied  flow  over  a  vertical  fence  when  the  motion  is  impulsively  started 
from  rest.  Solutions  have  been  calculated  for  the  following  two  cases: 

R  =  397,  P  =  10,  R.  =  1.58;  and  R  =  5000,  P  =  1,  R.  =  1.58.  The 
c  r  i  c  r  i 

results  are  presented  as  contour  plots  of  the  flow  variables,  generated 
directly  by  the  computer.  A  new  iterative  method  was  used  for  numerical 
integration  of  the  governing  equations  and  was  found  to  be  very  effec¬ 
tive.  The  method  is  based  on  the  combination  of  the  Crank-Nicolson 
method  with  a  strongly  implicit  iterative  method  developed  by  Stone  (1968). 
Full  details  of  the  method  are  given  together  with  a  summary  of  its  per¬ 
formance  on  the  problems  calculated.  It  is  shown  that  proper  treatment 
of  the  inflow  and  outflow  boundary  conditions  is  crucial  for  successful 
numerical  modeling  of  stratified  flows. 
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Chapter  I 
INTRODUCTION 

When  stably  stratified  air  flows  over  an  obstacle,  air  currents 
ter.i  to  fellow  a  wave-like  pattern  in  the  lee  side  of  the  obstacle. 

This  phenomenon,  called  mountain  lee-waves,  has  been  observed  in  the 
atmosphere  (Holmboe  and  Klieforth,  1957).  Lee-waves  have  beer,  studied 
analytically  with  the  help  of  perturbation  theory  for  steady  lee-waves 
of  small  amplitudes  (Lyra,  1943;  Queney,  1947;  Scorer,  1949).  For  a 
particular  boundary  condition  upstream,  /p  u  =  constant,  where  p  is 
the  fluid  density  and  u  is  the  streamwise  velocity,  a  nonlinear  model 
for  steady  lee-waves  of  finite  amplitude  was  investigated  by  Long  (1953), 
Yih  (1960),  Pao  (1967),  and  Miles  (1968,  1969). 

These  analytical  approaches  have  contributed  significantly  to  the 
understanding  of  lee-wave  phenomena.  However,  in  nature,  lee-wave  flows 
usually  involve  waves  of  finite  amplitude  and,  in  addition,  the  boundary 
conditions  are  complicated.  Moreover,  the  flow  may  be  in  a  transient 
state,  uence,  for  a  bettor  understanding  of  lee-wave  phenomena,  laboratory 
simulation  and  numerical  modeling  with  a  computer  are  the  only  approaches 
currently  available  apart  from  prototype  study  by  field  investigation. 
Laboratory  simulations  have  been  conducted  by  Long  (1959)  with  a  liquid 
model*  and  by  Lin  and  Binder  (1967)  in  a  wind  tunnel. 

In  laboratory  simulation,  one  usually  encounters  limitations  in  the 
variability  of  flow  parameters  such  as  flow  velocity  and  stratification. 

In  particular,  Reynolds  numbers  are  usually  quite  low  when  high  Froudc 
numbers  i«r  low  Richardson  numbers)  are  required.  Consequently,  molec¬ 
ular  effects  arc  inevitable  and  greatly  influence  the  flow  patterns  and 
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simulation  analysis  becomes  difficult.  As  far  as  numerical  computation 
is  concerned,  a  numerical  simulation  by  computer  is  usually  possible 
provided  a  stable  and  rapidly  convergent  scheme  is  available.  Neverthe¬ 
less,  a  numerical  computation  itself  is  by  no  means  a  trivial  and  straight¬ 
forward  manipulation.  In  particular,  when  a  nonlinear  partial  differen¬ 
tial  equation  such  as  the  Navier-Stokes  equation  is  involved,  a  numerical 
experiment  judged  with  physical  understanding  is  always  neeoed  because 
rigorous  analyses  of  computational  stability,  convergence,  and  error 
propagation  are  not  available  yet,  although  some  approximate  analyses 
based  on  the  existing  mathematical  theories  have  been  explored.  There¬ 
fore,  our  understanding  will  be  greatly  dependent  on  the  numerical 
experiment  coupled  with  these  approximate  analyses.  However,  even  though 
the  results  of  a  numerical  simulation  may  be  comparable  to  those  of 
laboratory  experiments,  a  numerical  simulation  is  still  desirable  since 
the  flow  parameters  can  be  varied  with  such  ease  and  then  more  information 
is  available.  Of  course,  one  must  recognize  that,  because  of  the  limita¬ 
tion  in  the  available  computer  storage  and  because  of  the  difficulty  in 
generating  turbulence  by  computer,  three-dimensional  problems  associated 
with  turbulence  must  still  be  studied  in  the  laboratory  or  by  prototype 
investigations. 

Numerical  studies  of  Ice-waves  have  been  carried  out  by  Hovormale 
C 1965 J  and  by  foldvik  and  Wurtelo  (1%?).  both  studies  treated  the 
transient  flow  of  inviscid  stratified  fluut  over  a  harrier  as  an  initial 
value  problem,  for  the  lack  of  suitable  boundary  conditions,  both 
upstream  and  downstream,  cyclic  boundary  condit:ons  were  ur.ed,  i.e., 
the  calculated  flow  conditions  downstream  of  the  barrier  were  fed  back 
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in  upstream  as  the  oncoming  flow  towards  the  barrier.  Some  important 
features  concerning  the  nonlinear  effects  on  lee-waves  were  obtained. 

In  the  numerical  experiment  described  in  this  report,  two-dimensional 
laminar  motion  of  an  incompressible,  viscous,  diffusive  and  nonhomogeneous 
fluid  over  an  obstacle  will  be  studied.  The  fruitful  results  from  the 
studies  by  liovermalc  and  by  Foldvik  and  Wurtele  indicate  that  treating 
a  transient  flow  instead  of  a  steady  flow  is  a  helpful  approach  from 
the  numerical  point  of  view  and,  thus,  in  this  numerical  study  a  transient 
flow  is  considered.  However,  we  shall  abandon  the  cyclic  boundary  condi¬ 
tion  and  make  a  detailed  study  of  the  effects  of  boundary  conditions  on 
the  computation.  In  addition,  a  new  iterative  method  is  introduced  in 
order  to  increase  the  speed  of  the  computation. 
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Chapter  II 

GOVERNING  EQUATIONS  AND  BOUNDARY  CONDITIONS 


2.1  Governing  Equations 

Consider  two-dimensional  laminar  motion  of  an  incompressible, 
viscous,  diffusive  and  nonhomoger.oous  fluid  in  a  uniform  gravitational 
field.  The  fluid  properties  other  than  the  density,  such  as  the  vis¬ 
cosity  viQ  and  thermal  conductivity  kQ  ,  are  assumed  to  be  uniform 
in  the  flow  field.  By  using  Boussinesq  approximation,  the  equations  of 
motion,  the  continuity  equation,  and  the  equation  of  energy  are  formu¬ 
lated  as  follows: 


3u  -  3u  -  3ii  3p 

p  — —  +  >i  —  +  w  — —  =  -  — + 

0  3t  3x  3z  3x 


MoV^U  , 


(2.1) 


-  ^  -  Pg  ♦  UqV2W  , 

3  z 


(2.2) 


3u/3x  +  3w/3z  =  0  ,  (2,3) 

and 

pc(^+ii^+w2i)=kV2i  ,  (2.4) 

°PUt  ax  3Z  0 

where  u  and  w  are  the  velocity  components  in  the  Cartesian  coordi¬ 
nates  x  and  z  ,  o  is  the  density,  po  is  the  reference  density, 
p  is  the  pressure,  C  is  the  specific  heat  capacity  at  constant 

r 

pressure,  T  is  the  temperature,  t  is  the  time,  and  v2  ■  32/3x2  ♦ 
32/3z2  is  the  Laplacian  operator.  In  Equation  (2.4)  the  heat  dissipa¬ 
tion  due  to  molecular  viscosity  is  neglected. 
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By  cross  differentiating  Equations  (2.1)  and  (2.2)  and  then 
eliminating  the  pressure,  one  arrives  at  the  vorticity  transport 
equation, 

^+u~+w~M==r-  %  g  +  vV2C  ,  (2.5) 
at  3x  3 z  /  o  3x 


where  v  =  mo/po  is  the  kinematic  viscosity,  the  vorticity  ?  is 
defined  by 

l  =  3w/ax  -  au/3z  ,  (2.6)* 

1  8T  1  3o 

and  =-  —  - - &•  has  resulted  from  use  of  the  equation  of  state 

*o  8x  po  3x 

for  the  incompressible  fluid  and  of  the  Boussinesq  approximation.  A 
stream  function,  i p  ,  which  satisfies  the  continuity  equation  (2.3)  is 
introduced  as 

u  =  3 ip/ 3 z  ,  and  w  =  -  3di/3x  ,  (2.7) 


and  the  vorticity  can  be  redefined  as 


Q  = 


V 1  iJj 


(2.8) 


For  convenience,  Equations  (2.5),  (2.7),  (2.8)  together  with  Equation 
(2.4),  presented  in  the  form 


3T 

-  3T 

-  37 

♦  u  — 

+  w  — 

at 

;.x 

,)Z 

(2.9) 


The  sign  convention  for  c  is  the  opposite  of  that  more  common ly 
used.  It  is  used  here  for  convenience. 
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in  which  k  =  ^/pC^  is  the  thermal  diffusivity,  will  be  normalized  to 
give  a  set  of  nondiuensional  equations, 
bet 

u  =  u/U  ,  w  =  w/U  ,  t  =  t/(L/U)  ,  ?  =  £/ (U/L)  , 

T-T 

x  =  x/L  ,  z  =  z/L  ,  T  =  4  =  uL/v  ,  Pr  =  v/k  , 

>p  =  iJi/UL  ,  and  F  =  U  =  ~ —  .  (2.10) 

/  ,  AT  vC 

\/gL  sr-  i 

*o 


U  ,  L  ,  and  AT  are  the  characteristic  velocity,  length,  and  temperature 
difference,  respectively,  and  Re  ,  Pr  ,  Fr  and  are  the  Reynolds 
number,  Prandtl  number,  Froude  number  and  Richardson  number.  The 
governing  equations  in  dimensionless  form  become: 


3t 


R.  |X 
1  3x 


(2.11) 


3T 

at 


+ 


+ 


P  R 
r  e 


72T 


(2.12) 


and 


C  8  -  V2ij/  , 


(2.13) 


where 

v2 


2 . 2  Governing  Equations  in  Finite  Difference  Form 

For  numerical  computation,  Equations  (2.11),  (2.12),  and  (2.13) 
will  be  discretized  into  stable  and  convergent  finite  difference  forms 
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and  then,  provided  the  initial  conditions,  boundary  conditions  and 
parameters  are  adequately  defined,  one  may  simulate  a  flow  with  the 
computer. 

In  Equation  (2.13)  the  derivatives  are  replaced  by  central  differ¬ 
ences  to  give, 


■  2<J/.  .  + 

_i lid _ 


LhJ. 


2i|». 


Jzis 


Ax 


Az' 


-  x..  •  (2.14) 

i.J 


where  Ax  and  Az  are  the  mesh  intervals  in  the  horizontal  and  vertical 
directions,  respectively,  the  subscripts  i  and  j  indicate  the  space 
coordinates  by  x  =  iAx  and  z  =  jAz  .  Figure  1  displays  the  arrange¬ 
ment  of  the  mesh.  It  can  readily  be  shown  that  in  respect  to  Equation 

(2.13) ,  Equation  (2.14)  has  a  truncation  error  e  »  0(Ax2)  *  0(Az2) 

where  the  notation  0  denotes  that  lim  ^7^-  is  finite.  Thus,  the 

6-+o  6 

smaller  the  mesh  size  is,  the  smaller  the  truncation  error  will  be. 

Discretization  of  Equations  (2.11)  and  (2.12)  presents  greater 

difficulty  since  one  must  consider  questions  of  convergence  and  of 

stability.  If  a  finite  difference  scheme  has  the  property  of  stability, 

the  difference  between  the  solution  of  the  finite  difference  system, 

say  0.  .  ,  ,  and  the  exact  solution  of  the  differential  system,  say 
1 » J  >  * 

0(iAx,jAz,kAt) ,  will  remain  bounded  as  time  increases  indefinitely  for 
fixed  mesh  sizes,  Ax  ,  Az  and  time  increment  At  .  Convergence  guaran¬ 
tees  that,  for  t  fixed,  0.  .  .  will  converge  to  e(iAx, jAz.kAt)  as 

*  *  J  »  * 

Ax  ,  Az  ,  and  At  approach  to  zero.  A  further  difficulty  arises  from 
the  nonlinearity  involved  in  the  coupling  between  Equations  (2.11)  and 

(2.13) .  Experience  shows  that  the  nonlinearity  can  produce  instability 
in  the  computations  although  stability  analysis  based  on  the  locally 
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linearized  finite  difference  system  indicates  a  stable  scheme.  For  this 
reason  numerical  experiments  are  always  needed. 

Generally  speaking,  there  are  two  classes  of  schemes,  explicit  and 
implicit,  for  solving  a  transient  problem.  In  explicit  schemes  the 
solution  at  a  given  grid  point  for  the  new  time  level  can  be  calculated 
from  known  quantities  at  the  present  or  previous  time  level.  On  the 
other  hand,  implicit  schemes  involve  the  solution  of  a  set  of  simultane¬ 
ous  equations  for  calculating  the  solution  at  all  the  grid  points  at 
the  new  time  level.  Obviously,  explicit  schemes  involve  less  algebraic 
manipulations  than  do  implicit  schemes.  However,  for  a  linear  system, 
stability  analysis  shows  that  the  time  increment  At  for  explicit 
schemes  is  restricted  and  usually  At  is  kept  small  to  satisfy  require¬ 
ments  of  stability  and  of  convergence;  whereas,  in  implicit  schemes  the 
time  interval  has  no  restrictions  except  those  imposed  by  considerations 
of  accuracy.  Consequently,  a  larger  time  increment  can  be  used  with 
implicit  schemes  and  this  can  result  in  less  computer  time  being  required 
than  in  the  case  of  explicit  schemes.  In  particular,  if  we  are  chiefly 
interested  in  the  asymptotic  solution  of  a  transient  flow,  implicit 
schemes  appear  to  be  preferable  to  explicit  schemes.  Here  we  would 
point  out  that  the  computational  scheme  which  is  best  for  i  given  prob¬ 
lem  will  often  be  determined  by  the  nature  of  the  physical  phenomena. 

For  instance,  if  a  flow  involves  phenomena  of  short  period  or  high 
frequency  the  time  increment  must  be  considerably  shorter  than  this 
period.  An  example  is  the  numerical  calculation  of  vortex  shedding 
phenomena  behind  an  obstacle  (Fromm  and  liarlow,  1963;  Thoman  and 
Szewczyk,  1966).  And,  of  course,  the  mesh  size  must  be  considerably 
less  than  the  characteristic  length  of  the  shedding  vortex,  ,'s  for 
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the  present  numerical  investigation,  we  shall  use  an  implicit  scheme 
since  transient  lee-waves  propagate  slowly  and  have  a  characteristic 
length  comparable  to  or  larger  than  the  obstacle  height  or  width. 

Even  when  an  implicit  scheme  is  used,  the  nonlinearity  of  the 
equations  can  cause  difficulties  when  the  Reynolds  number  Re  becomes 
large.  If  the  convective  terms  of  Equation  (2.11)  are  discretized  with 
central  space  differences,  the  off-diagonal  terms  of  the  amplification 
matrix  of  the  difference  equations  can  become  large  in  magnitude  com¬ 
pared  to  the  principal  diagonal  te’ms  as  Re  increases  and,  then,  the 

absolute  values  of  the  eigenvalues  of  tne  amplification  matrix  may  be 

greater  than  unity.  This  situation  has  been  improved  somewhat  by  use 
of  the  alternating  direction  implicit  (ADI)  method  (Peaceman  and 
Rachford,  1955)  for  moderate  Reynolds  numbers  (Pearson,  1965).  For 
higher  Reynolds  numbers,  the  backward  and  forward  difference  method 
for  convective  terms  has  proved  to  be  very  effective  (Thoman  and  Szewczyk, 
1966).  This  method  was  first  suggested  by  Lelevier  (Richtmyer  and 
Morton,  1967)  and  consists  in  representing  d^/Zx  (or  3$/3z)  by  the 

backward  or  forward  finite  difference  depending  on  whether  u  (or  w) 

is  positive  or  negative.  In  this  way,  the  amplification  matrix  is  made 
positive  definite  and  diagonally  dominant;  thus,  the  absolute  values  of 

the  eigenvalues  of  the  amplification  matrix  are  less  than  unity  and  the 
stability  criterion  becomes  weakly  dependent  on  the  Reynolds  number. 

This  method  together  with  the  ADI  method  has  been  used  by  Pao  (1969)  for 
calculation  of  homogeneous  flow  of  a  viscous  incompressible  fluid  over  a 
horizontal  flat  plate  of  finite  length. 

In  the  present  study,  for  high  Reynolds  number  this  backward  and 
forward  diffeicncc  method  is  applied  to  the  convective  terns  and  then 
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the  Crank-Nicolson  implicit  scheme  is  used.  The  governing  equations 
in  finite  difference  form  are  formulated  as  follows: 

(i)  Temperature  field: 


(ii)  Vorticity  field: 


12 


and 

(iii)  Stream  function  field: 


2n, 

,  + 

.  , 

i-l.J 

I 

+  ■ 

Ax2 

= 

9 

in  which 

(a) 

A2  * 

1  , 

A3 

=  0 

for 

(b) 

A2  = 

o  , 

A3 

=  1 

for 

(c) 

A4 

1  , 

As 

=  0 

for 

and 

(d) 

A4 

0  , 

A5 

=  1 

for 

In  Equations  (2 

!.  15) 

and 

(2.16) 

,  the 

the 

time 

steps 

(n+1) At 

and 

nAt  . 

ib.  .  ,  -  2\b.  .  +  \p.  . 

Ti  1^1  " 1  1  T  1  1  . 


A  Z 


for  it?  .  >  0  , 


n 

i.  . 

i»J 

n 


(2.17) 


(2.18) 


quantities  at  the  nth  time  step  are  known,  while  those  at  the  (n+1) 
time  step  are  unknown  and  must  be  calculated. 


2.3  Initial  and  Boundary  Conditions 

In  the  atmosphere,  flow  conditions  may  vary  from  time  to  time  and 
we  may  use  observed  data  as  the  initial  conditions  and  calculate  the 
flow  as  it  develops  from  that  state.  However,  this  makes  the  problem 
complicated  and,  in  this  numerical  experiment  we  consider  a  simpler 
condition.  Consider  a  vertical  fence  of  infinitesimal  thickness  and 
height  2L  in  a  fluid  of  unlimited  extent.  The  fluid,  initially  uni¬ 
formly  stratified  in  the  vertical  direction  and  at  rest,  is  impulsively 
accelerated  ii:  the  horizontal  direction  so  that  the  velocity  far  from 
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the  fence  is  U  and,  thereafter,  the  velocity  is  maintained  constant. 

The  flow  is  initially  irrotational  and  the  vorticity  which  is  generated 
as  a  result  of  the  no-slip  boundary  condition  on  the  vertical  fence  is 
concentrated  only  on  the  fence.  As  time  passes,  vorticity  spreads  into 
the  flow  field  under  the  combined  action  of  viscous  diffusion  and  of 
convection.  The  temperature  distribution  at  the  initial  instant  when 
the  flow  starts  impulsively  is  linearly  distributed  with  height  and  the 
temperature  is  constant  along  a  streamline;  the  characteristic  tempera¬ 
ture  difference  is  defined  by  AT  =  T^  -  Tq  where  T^  is  the  tempera¬ 
ture  far  upstream  at  height  L  measured  from  the  center  of  the  vertical 
fence.  As  time  passes,  the  temperature  distribution  will  be  changed 
as  a  result  of  thermal  diffusion  and  of  convection. 

As  far  as  the  boundary  conditions  are  concerned,  only  the  no-slip 
condition  for  a  viscous  fluid  is  clearly  defined  on  the  boundary  of 
obstacle.  However,  the  governing  equations  are  expressed  in  terms  of 
vorticity  and  the  velocity  is  in  turn  calculated  from  the  stream  function 
equation;  in  other  words,  there  is  no  a  priori  boundary  condition 
available  for  vorticity.  In  the  following,  we  express  the  vorticity 
on  the  obstacle  boundary  in  terms  of  the  vorticity  and  stream  function 
in  the  neighborhood. 

Referring  t.)  figure  2,  we  express  the  stream  function  at  A 

in  terms  of  the  stream  function  ^  and  its  derivatives  at  B  by  a 
Taylor  series  expansion 


%  ♦  Ax  S 


/  .)ni 


>  JX 


Ax' 


I  ]  ♦  Ax3 
I  ?  31 

B 


<)x 


aVi 

i)x4 '  B 


♦  0(Ax'J) 


(2.19) 
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Equation  (2.19)  can  be  related  to  the  vorticity  and  its  derivatives  at 
B  by  virtue  of  the  no-slip  boundary  condition  at  B  as 


Ax2 


f&j  . 

Ax4  / 32C 

32;  I 

Mb 

4!  lax2' 

az2l 

which  is  then  rewritten  in  the  form 


‘B  ■  77  (VV  -  k  SA  *  TT  <^>8  * 

A.X 


'  77  -  7  'a  *  TT  Re  lli)B  *  °(“3’  ‘2'20> 


by  using  the  Taylor  expansion  of  about  B  and  the  no-slip  condition. 
The  same  derivation  was  given  by  Hung  (1966).  Now,  if  (3s/at)g  is 
expressed  by 


n+1 


-  Cr 


At 


and  the  remaining  terms  of  the  right  hand  side  of  Equation  (2.20)  are 
given  the  values  corresponding  to  the  nth  time  step,  then  we  have 


n+1 


8At 

R  Ax5 
e 


3  ,n-.  _  1  n  .  „n 

”  2  2  CA  C3 

Ax 


1  + 


R  AX2 

0 

8At 


+  O(AxAt)  +  0(At2)  . 


(2.21) 


For  the  vorticity  on  the  upstream  face  of  the  vertical  fence,  a  similar 
equation  is  available: 


IS 


n+i 

V 


8At 

R  Ax2 
e 


3  ,.n  ,n  .  1  n  n 

2  ^B'  "  2  ^A’  C 

Ax 


B' 


1  + 


R  AX2 
e 


8At 


♦  O(AxAt)  +  0(At2)  . 


(2.22) 


Hence  the  vorticity  on  the  boundary  of  the  vertical  fence  at  the 
(n+1)  time  step  can  be  evaluated  from  the  quantities  at.  the  nth  time 
step,  and  then  is  used  as  a  boundary  condition  for  the  calculation  of 
vorticity  at  interior  points  at  the  n+1  time  step.  In  the  case  when 
an  iterative  process  is  used  for  the  calculation  of  vorticity  at.  a 
given  time  step,  and  at  the  (m+1)  iteration  are 

calculated  from  the  values  of  the  mth  iteration  by. 


.(m+1) 


(*■ 


(m) 


Ax 


'  *Am))  - 


(2.23) 


and 


r  (m+l) 
’B' 


Ax 


Up  -  *P>  ■  3  4?  *  °(“2) 


(2.24) 


Subsequently  they  arc  used  as  boundary  conditions  for  calculating  vor¬ 
ticity  at  interior  points  at  the  (m+1)  iteration. 

As  for  the  boundary  value  of  vorticity  on  the  comer  of  the  verti¬ 
cal  fence,  C  ,  a  different  expression  is  derived  as  follows.  Expanding 
the  vorticitics  c.j  ,  r>l;,  ,  c0  »  and  in  Taylor  series  about  C  , 
we  have 

■2tu  *  Cf  •  -  CC  ♦  (^f)c  •  0(44’)  . 

and 


l 
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As  far  as  temperature  distribution  on  the  obstacle  is  concerned,  we 

consider  a  uniform  distribution  for  convenience,  i.e.,  f  =  T  or 

o 

T  a  0  . 

On  the  outer  boundary  beyond  ;he  obstacle  the  boundary  conditions 
are  not  defined  a  priori.  Upstream  very  far  from  the  obstacle,  one  may 
define  a  prescribed  boundary  condition.  However,  because  of  the  pres¬ 
ence  of  disturbances  propagating  from  the  obstacle  both  upstream  and 
downstream,  the  boundary  conditions  for  both  inflow  and  outflow  bound¬ 
aries  have  to  be  defined  in  a  less  restrictive  way  by  an  extrapolation 
method  derived  from  the  Milne  predictor  formula. 
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yi+1  =  yi-3  +  T  (2yi  "  yi-i  +  +  °^x5)  »  C2-28) 

where  the  prime  denotes  the  derivative  with  respect  to  x  and  y  is 
a  function  of  x  .  By  using  the  formulae 

n  =  i h  -  16),i-3  *  36>-i.2  -  *  25yi>  *  • 

i'l-i =  *  6^i-3  -  *  ^i-i  *  3>V  *  °(4x'‘)  • 

and 

yi-2  =  lh  (yi-4  '  8yi-3  *  Syi-1  -  yi>  *  • 

Equation  (2.28)  is  formulated  as 

yi+l  =  yi-4  "  Syi-3  +  10yi-2  -  10yi-l  +  5yi  +  0(Ax5)  ‘  (2*29^ 

Hence,  with  an  error  of  order  Ax5  ,  the  stream  function,  vorticity, 
and  temperature  on  the  outflow  boundary  are  expressed  by  extrapolating 
from  the  neighboring  upstream  points  as  follows: 


h.l.j  -  *i-4,j  -  S*i-3,j  *  10*i-2.j  -  10h-lJ  *  5*i,j  -  ‘2-30a> 


h.l.j  ■  4i-4„i  -  5<i.3,j  *  10'i-2,j  -  *  5S.j  •  (2-30W 


and 


*  Ti-4, j  -  5Ti.3.j  *  “WJ  - 


18 


A  similar  but  different  form  was  used  by  Hung  (1966)  successfully  for 
a  laminar  flow  in  a  conduit  expansion.  Similarly,  these  variables  on 
the  inflow  boundary  are  expressed  by  extrapolating  from  the  neighboring 
downstream  points 


♦i-l.i  =  *  10+i.2,j  '  “VlJ  *  5*i,j 


”  Ci*4,j  "  5ci*3,j  *  10ti*2,j  ■  *  5ci,j 


and 


(2.31a) 

(2.31b) 


(2.31c) 


The  detailed  study  of  boundary  conditions  Equations  (2.30)  and  (2.31) 
will  be  presented  when  the  computational  results  are  described. 

For  an  exact  simulation  of  flow  over  an  obstacle  by  computer,  the 
computation  should  cover  the  whole  flow  field  1234  (see  Fig.  3),  espe¬ 
cially  when  some  asymmetric  flow  patterns  such  as  shedding  vortices 

Ri^ 

are  present.  These  would  be  expected  to  occur  when  o  =  . . —  <  1  , 

log  Re/40 

and  Rg  >  40  (Pao,  1968).  However,  because  of  the  limitation  of  com¬ 
puter  capacity,  a  flow  field  of  half  space  1256  was  considered  instead, 
it  is  assumed  that  the  flow  field  is  symmetric  about  the  line  56. 
Actually  this  is  not  a  bad  assumption  when  o  >  1  according  to  Pao. 
Hence,  we  assume  that  when  t  >  0 


♦l.l  *  0 


'i , 1 


■  0 


T.  ,  =  0 

1,1 


or 


T 

o 


(2.32) 


for  the  boundary  conditions  on  the  lino  of  symmetry,  5b,  and 
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wi,N  =  0  * 


9/9x(T._n)  =  0 


(2.33) 


for  the  upper  boundary  condition  on  the  line  12 
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Chapter  III 
ITERATION  METHOD 

In  the  numerical  solution  of  Equations  (2.15),  (2.16),  and  (2.17) 
an  iterative  process  is  always  encountered  because  of  the  boundary  con¬ 
ditions  and  the  nonlinearity  of  the  equations.  Iterative  methods,  such 
as  the  Gauss-Seidel  process,  SOR  (successive  over  relaxation) ,  or  ADI 
(alternating  direction  implicit  method)  may  be  used.  However,  in  this 
study  a  new  strongly  implicit  iterative  method  proposed  by  Stone  (1968) 
will  be  used.  This  method  has  several  advantages  over  the  methods 
previously  used.  Firstly,  its  rate  of  convergence  does  not  depend 
strongly  on  the  nature  of  the  coefficient  matrix  of  the  equations  to 
be  solved;  secondly,  it  is  not  sensitive  to  the  choice  of  iteration 
parameter;  and  thirdly,  it  reduces  significantly  the  computational 
effort. 

3. 1  Statement  of  Problem 

There  are  many  methods  available  for  solving  a  set  of  linear 
algebraic  equations,  expressed  in  matrix  form  as, 

[M]f  *  q  .  (3.1) 

Direct  elimination  is  always  possible  but  it  is  time  consuming  and 
requires  excessive  computer  core  if  the  dimension  of  the  matrix  is 
large.  Methods  based  on  triangular  resolution  of  [M]  into  a  lower 
and  an  upper  triangular  matrix,  [l]  and  [IJj  respectively,  can  result 
in  significant  savings  in  computational  effort  and  in  storage  require¬ 
ments  as  compared  to  other  direct  methods  of  solution  when  [Mj  is  a 
band  matrix.  In  the  particular  case  when  (MJ  is  tri-diagonal  the 
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savings  are  very  gicat  indeed.  When  partial  differential  equations  in 
two-dimensional  space  are  approximated  by  implicit  finite  difference 
equations  involving  function  values  only  at  the  reference  mesh  point  and 

at  its  nearest  four  neighbor  points,  the  matrix  of  coefficients  of  the 

difference  equations  has  non-zero  elements  only  along  five  diagonals, 
viz.,  the  principal  diagonal  and  that  on  either  side  of  it  and  two  off- 

diagonal  lines.  In  general,  when  such  a  matrix  is  resolved  into  [L] 

and  [U]  ,  these  matrices  do  not  preserve  the  sparse  structure  of  [M] . 
However,  if  a  matrix  which  approximates  [M]  can  be  found  such  that  it 
is  the  product  of  a  lower  and  of  an  upper  triangular  matrix,  each  of 
which  has  non-zero  elements  on  only  three  diagonal  lines  corresponding 
to  the  non-zero  diagonal  lines  in  [M]  ,  then  it  is  possible  to  solve 
the  set  of  equations  by  an  iterative  method  and  to  preserve  the  benefits 
arising  from  the  sparseness  of  the  matrix  [M] . 

The  expressions  presented  in  Equations  (2.15),  (2.16),  and  (2.17), 
can  be  simplified  in  a  general  form 


B.  .  T.  .  , 
i»J  i.J-1 


+  l). 


T. 


i»j  i-1 > j 


+  E.  .  T.  .  +  F.  .  T.  . 
i.J  i.J  i,J  i*l. J 


*  H 


i»j  i,j+l 


*  q 


i,3 


(3.2) 


which,  expressed  in  matrix  notation,  has  the  form 
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(3.3) 

or 

[M]T  »  q  ,  (3.1) 

where  the  matrix  [M]  is  a  (I+l)x(J  +  l)  square  array.  The  elements  of 
the  matrix  are  composed  of  the  coefficients  of  temperature,  vorticity, 
and  stream  function  in  the  respective  Equations  (2.15),  (2.16),  and 
(2.17).  T  is  a  vector  composed  of  the  sequenced  temperatures,  or  vor- 
ticities,  or  stream  functions  over  the  computation  field,  q  is  a 
vector  composed  of  the  sequenced  known  quantities.  In  the  matrix  [M] 
non-zero  elements  occur  only  along  the  five  diagonal  lines  shown. 


3.2  Algorithm  for  Iteration 


In  order  to  solve  Equation  (3.2)  iteratively,  an  inductive  algorithm 

instead  of  the  algorithm  described  by  Stone  (1968)  will  be  introduced. 

In  principle,  both  algorithms  have  the  same  outcome,  but  the  former,  it 

seems  to  us,  has  a  clearer  and  simpler  derivation. 

Assume  that  T.  .  can  be  evaluated  from  T.  .  ,  and  T.  .  , 
i.J  i,J+l  i+l.J 

which  are  already  known.  Then  let 


T.  .  =  e.  .  T.  ,  .  ♦  f.  .  T. 
i.J  i.J  1+1»J  i.J  1 


in  which  e.  .  ,  f.  .  ,  a.id  d.  . 

i.l  i.J  i.J 

by  Equation  (3.4),  ^  ^  and  T 

pressed  in  the  forms 


.  ,  +  d.  . 
»J  +  1  i.J 


(3.4) 


are  subject  to  determination.  Hence 

.  .  .  which  are  unknown  can  be  ex- 
i.J-1 


T. 


i-l.j 


■  e.  .  .  T.  . 
i-l. J  i.J 


fi-l,j  Ti-l,j+l  +  di-l,j 


(3.5) 


and 


i.j-l  Ci.j-1  'i+l.j-l  ’  Ai,j-1 


T 


♦  f. 


T.  .  ♦  d 
i.J 


i.j-l 


(3.6) 


Substitute  Equations  (3.5)  and  (3.6)  into  Equation  (3.2),  and  then  we 
have 


T. 

i.J 


I-.  . 

_ Id _ _ 

D.  c.  ,  >E.  .♦B.  .{.  .  , 

i.J  i-l. J  i.J  i.J  i.J-1 


T. 


i  +  l.  j 


“i.  t.J  ‘-i.n-i.j.i 


B.  .e.  .  ,T.  .  .  , 
i,j  i.j-l  i*l,; i-l 


(3.7) 


Equating  the  right  hand  sides  of  Equations  (3.4)  and  (3.71,  we  obtain 
the  following  relations: 
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F. 

_ LJ _ _ _ 

D\  .e.  ,  .  +  E. .  +  B.  7f.  TT 
i.J  i-l, J  i,2  1.3  1,3-1 


=  e 


i,j  ’ 


(3.S) 


H.  . 

- Lij _ =  f 

I).  .e.  ,  .  +  E.  .  +  B.  .f.  .  .  i ,  j  * 
i,3  1-1,3  i,3  i,3  1,3-1 


(3.9) 


.  -  D 


-  B.  .d.  . 


D.  .e.  ,  . 
i,3  1-1,3 


,  ♦  li,  .  +  B.  .f.  .  .  Qi, j  ' 


i,3 


i,3  i.3-1 


(3.10) 


and 


U.  .f.  ,  .  =0  , 

1.3  i-l,3 


B.  .e.  .  .  =  0  . 

i,3  1,3-1 


(3.11) 


(3.12) 


In  the  expressions,  Equation  (3.8)  through  Equation  (3.12),  the  capital¬ 
ized  elements  are  all  known.  If  the  elements  with  subscripts  (i-l,j) 
or  (i,j-l)  are  also  known  from  previous  calculations,  then  in  the  above 

five  relationships  we  have  only  three  unknowns,  e,  .  ,  f.  .  ,  and  d.  .. 

i,3  i,3  i,3 

In  other  words,  these  three  unknowns  are  indeterminate  if  all  five 

relationships  have  to  be  satisfied;  equivalently,  we  can  say  that  there 

is  no  way  to  resolve  the  matrix  [M]  of  Equation  (3.1)  exactly  into  a 

product  of  lower  and  upper  triangular  matrices  with  the  same  sparse 

structure  as  [it]  ,  as  we  previously  noted. 

In  Equation  (3.7),  the  terms  containing  T\  ^  and  ^  j 

are  the  source  of  the  trouble.  However,  if  T.  ,  .  ,  and  T.  .  .  , 

x-l,j+l  i+l, 3-1 

can  be  expressed  approximately  in  terms  cf  T\+^  ^  ,  or  T.  ^  ,  or 

even  T.  .  and  other  known  quantities  when  we  evaluate  T.  .  ,  then 
*  *  J  1  »  J 

we  obtain  three  relationships  for  three  unknowns.  Expanding  T.  ,  .  ,  , 

1*1  p  Ji  *  I 

T.  j  .  ,  and  T.  ^  in  Taylor  series  about  the  grid  point  (i,j),  we 
have 
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T. 

i-l.j+1 


+  T 


i.j  +  l 


+  T.  .  .  +  O(AxAz) 


(3.13) 


Similarly, 


(3.14) 


Now  it  is  realized  that  the  right  hand  sides  of  the  expressions  (3.13) 
and  (3.14)  have  a  truncation  error  O(AxAz)  and  thus  an  iterative  process 
is  necessary  for  the  solution  of  Equation  (3.3).  Therefore,  we  introduce 
an  iteration  parameter  a  for  both  expressions  (3.13)  and  (3.14)  in  the 
forms 


T.  .  .  .  ■  ot(-T.  .  +  T.  .  ,  +  T.  ,  .) 
i-l.  1  +  1  i.J  i.J  +  1  i-l, ] 


(3.15) 


and 


T. 


i+l.j-1 


of-T. 


i.J 


+  T.  ,  .  +  T.  .  .) 

l+l, j  i,j-l 


(3.16) 


In  principle,  two  iteration  parameters,  say  a  and  6  may  be  introduced 
to  the  expression  Equations  (3.13)  and  (3.14)  respectively.  However,  we 
use  only  one  iteration  parameter  a  for  convenitn.ee  and  for  simplifying 
the  numerical  experiment. 

First  consider  the  expression  Equation  (3.15).  By  virtue  of  the 
Equation  (3.5),  Equation  (3.15)  becomes 


•l.j  1 i-l  ,j*l  +  Ui-l,j^ 


» 


which  gives 
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(-1  +  e.  .  .)T.  .  +  T.  .  .  +  d.  .  . 

T.  ,  .  .  =  a  - bhl-bO - bill - i ihl 

i-l, 3+1 


(3.17) 


l-o  f 


i-l,j 


Similarly,  by  use  of  Equation  (3.6),  T.  .  .  .  is  expressed  in  the 

l+i, j-i 

form 

(-1  +  f.  .  ,)T.  .  +  T.  .  .  +  d.  . 


T.  ,  .  ,  =  a 

l+l, j-1 


i,j-l  i.j  i»l.j  i, 3-1 


1  -  a  e. 


(3.18) 


i » j-1 


Finally,  substituting  Equations  (3.17)  and  (3.18)  into  Equation  (3.7) 
the  result  is  obtained, 
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i,J 


i.J 
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-bibllLtl 
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1-af 
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1-ae.  .  , 
i.J-1 
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i,j  l-«ejj.i 
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D . . (e. ,  ,-af.  .  .)  B.  ,(f.  .  ,-ae.  .  ,)  i+l.j 
h  +  ~bJ — LlijJ - i"  *„»  J.  +  i,J  i,j~l  i.J-1 


qi.j  ■ 

■  "i.j  > 

di-l.j 

I  B.  .  ) 

l.J 

di.j-l 

|l:reLi-i  1 

u.  (e.  ,  .-t,f.  ,  .) 

t  ♦  x>y  blxi  *-*>j  . 

I.J 


i-l.j 


.j-1 


(3.19) 


By  comparing  Equations  (3.4)  and  (3.19),  e.  .  ,  f.  .  ,  and  d.  are 

!*J  i.J  i.j 

obtained  inductively  as  the  following  expressions. 


r-r-K' 
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aB.  .e.  .  . 

F  + 
i.j  1-012, 


6i»j 


E. 


i.J 


U.  .(e.  .  .-af.  ,  .)  §"!  .  (f .  .  ,-ae.  .  .) 

+  X»J  1"1»3  +  _ L»J..lA _ 1,J  iL- 


1-af 


i-1  ,3 


1-ae. 


i,j-l 


(3.20) 


aD.  .f.  . 

1lL  Izhl 


D.  .(e. 


H.  .  +  ■-■■■4 

i,i  1-af.  .  . 

J _ i-l.J 

s,  ,  ,  -af,  ,  1 


.)  B.  .(f. 

i  *  i  i  v  i 


E.  .  ♦ 

i.J 


jv  i-i,j  i-i, y  +  i,]v  i,j-i  i,j-i 


1-af 


i-l.j 


1-ae 


i.j-l 


(3.21) 


and 


d.  . 
i.J 


1  Dwi  |d  . 

1  B.  . 

1>J  d 

qi,j  '  | 

1-ae.  .  .  i,j-l 

1  1,3-1/ 

E .  .  ♦  ■ 
1.3 

D.  .(e.  .  .-af.  ,  •) 
1,3  i-l.J  1-1.3 

B.  .  (f.  .  ,-ae.  .  . 
i»3  1.3-1  1,3-1 

1-af.  .  - 
1-1.3 

‘-"i.j-i 

(3.22) 


Now  the  computational  procedure  is  straigntforward  and  is  described 
as  follows: 

(1)  Since  all  the  capitalized  elements  and  q.  .  are  known  quanti- 

1 ,3 

ties  at  all  grid  points,  aid  since  Equations  (3.20),  (3.21),  and  (3.22) 

give  the  expressions  of  e.  .  ,  f.  .  ,  and  d.  .  in  terns  of  these  at 

1 » J  *  >  J  1  >  J 

grid  points  (i,j-l)  and  (i-l,j),  one  can  simply  evaluate  e.  .  , 

1  >  J 

f.  .  anu  d.  .  in  the  ^holc  field  starting  from  the  grid  point  at  the 
left,  lower  corner  of  the  field  to  that  at  the  right,  upper  corner  of  the 
field. 

(2)  When  ail  the  c.  .  ,  f.  .  ,  and  d.  .  have  been  calculated, 

*  >  )  *  *  J  1  *  J 

T.  can  be  calculated  "explicitly"  by  Equation  (3.19)  from  the  right, 

*  *  J 

upper  comer  of  the  field  to  the  left,  lower  corner  provided  suitable 
boundary  conditions  arc  given. 
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It  is  interesting  to  note  here  that  the  original  set  of  Equations,  (3.2), 
is  an  implicit  one  and  that,  by  applying  the  above  stated  algorithm,  one 
can  solve  the  equation  (approximately)  in  an  "explicit"  manner. 

Since  the  right  hand  sides  of  Equations  (3.15)  and  (3.16)  are  only 
the  approximate  expressions  for  T\  ^  ^+(,  and  T\+^  ^  ^  with  trunca¬ 
tion  error  0(AxAzj ,  the  calculations  from  the  expressions  Equations 
(3.19)  to  (3.22)  give  only  a  first  approximation  to  the  exact  solution, 
and  an  iterative  process  is  needed  to  improve  the  accuracy.  Consider 
Equation  (3.1)  at  the  (m+1)  iteration. 


M  =  <5 


(3.23) 


which  can  be  rewritten 


[M]  (T(lllM)  -  TW)  =  q  -  [M]  T( 


(3.24) 


Let  -  Tl‘m^  ,  the  changes  of  T  from  iteration  m  to 

m+1,  and  =  q  -  [M]T^I;  ,  the  residuals  at  iteration  m  .  There- 

(ni) 

fore,  replacing  q,  .  in  Equations  (3.22)  and  (3.19)  by  R.  .  , 

*  i  )  1  I  J 


=  q.  .  -  (B.  ,T  . 
1.3  1*3  1.3- 


,T<n).  +  D.  .T^m]  .  +  E.  .T^  +  F.  .T?m} 

l  i,j-l  i,j  i-l.j  i,3  1,3  i»3  i+l,3 


+  li.  .  .)  , 

i,3  i,i+r 


(3.25) 


and,  replacing  T.  .  ,  T.  .  .  ,  and  T.  .  .  in  Equation  (3.19)  by 
1  *  J  1  >  J  *  *  *  *  A  i  J 

,  and  ,  respectively,  we  obtain  the  algorithm 

fo;  the  iterative  solution  ui'  Equation  (3.1): 


.  c.  .si™;1’  .  t.  .  «<“*'>  .  a.  . 

1,3  i.J  1+1,3  i,J  1,3+1  1.3 


(3.26) 
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where  e.  .  and  f.  .  have  the  same  definitions  as  in  Equations  (3.20) 
i.l  i,l  ” 

and  (3.21),  while  d.  .  is  now  given  by, 

i» 1 


d.  . 
i.l 


rP'P  - 
1,1 

D.  .  \ 

*,1  d 

B.  .  \ 

1,1  d 

1-af.  .  .  1  i-l, j 
i-l,l  I 

1-ae.  .  .  /  ai,j-l 
1,1-1  1  J 

P 

D.  .(e.  .  .-af.  .  .) 
x  ill  1-1,1 _ i-l,l 

B.  .(f.  .  ,-ae.  .  J 
i,l.  1,1-1  i,l-l' 

C  •  • 

1-af.  1  . 

1-1,1 

1-ae.  .  , 

1,1-1 

(3.27) 


The  iterative  process  can  be  continued  until  the  residuals  rP1'} 

satisfy  a  certain  criterion,  and  then  tP‘1+^  is  the  solution  at  the 

1 » 1 

corresponding  grid  points. 


3.3  Convergence  Rate 

The  solution  expressed  in  Equation  (3.19)  together  with  Equations 
(3.20),  (3.21),  and  (3.22)  is  an  approximate  solution  to  Equation  (3.1). 

[M]  T  =  q  ,  (3.1) 

but  is  the  exact  solution  to 


B.  .T.  .  ,  +  U.  .T.  .  • 
1,1  x, j-1  i,j  i-l, J 


+  E.  .T.  .  +  F.  .T.  .  .  +  H.  .T.  .  , 

i.J  i,l  1,1  1+1,1  i,l  1,1+1 


-  D.  .f.  .  .[T.  ,  . 
i,l  i-l, ll  i-l, 1  +  1 


a(-T.  .  +  T.  .  .  +  T.  ,  .)] 

1,1  1,1  +  1  i-l,  y1 


(5.28) 

The  proof  can  easily  be  obtained  by  substituting  Equations  (3.5)  and 
(3.6)  into  Equation  (3.2)  and  by  the  fact  that  the  resultant  equation 
is  identical  to  the  result  of  combining  Equations  (3.15)  and  (3.16)  with 
Equation  (3.7).  The  modified  Equation  (3.2)  in  matrix  form  is 
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|[M]  +  [N]:  T  =  q  .  (3.29) 

Therefore,  the  convergence  rate  for  the  corresponding  iterative  process 
for  Equation  (3.29),  i.e., 


j[M]  +  [N])  T(W+1)  =  | [Mj  +  [N])TC,n)  +  (q  -  [M]  T(m)j  (3.30) 

can  be  analyzed. 

If  T  is  the  exact  solution  to  Equation  (3.1),  then  the  error  of 
at  the  mth  iteration  is  defined  by 

|(m)  =  j(m)  .  j  4  (3.31) 

Thus,  Equation  (3.30)  can  be  rewritten  into  an  error  equation 
|[>1]  •  =  |[M]  .  [N]jJW  -  [M]t« 

.  [Nlt(m)  ,  (3.32) 


which  can  also  be  expressed  in  the  form 

p  (ln+D  +  n  e(»+1)  ♦  p  r(m+D  +  f  (m+  i)  +  u  _>+D 

r,  r  r  (m+1)  ,  (m+1)  (m+1)  (m+1),.-, 

_  ,  (m+1)  ,  (m+1)  (m+1)  (m+l)..-. 

"“i.j'i.j-l  ei*l,j-l  ' 

-  -D  f  re(m)  -  e(-e^  *  +  r!B>  .)] 

"  Ui,j  i-l,jl  i-l,j+l  1  i,j  i,j  +  l  i-l,  rJ 


„  r  (w)  ,  (m)  ^  „(m) 


(3.33) 


In  general,  the  coefficients  in  Equation  (3.33)  are  not  constants 
but  vary  from  one  grid  point  to  another,  from  one  time  step  to  the 
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next,  and  from  one  iteration  to  the  next.  However,  for  indicating  how 
the  convergence  rate  behaves  we  simply  assume 


B.  .  =  li. 

i,l 


B 


and 


E.  .  =  -  2(B  +  D) 

1,1 


(3.34) 


These  conditions  may  satisfy  the  stream  function  Equation  (2.17)  pro¬ 
vided  Ax  =  Az  .  Further,  we  assume 


-D.  .  f.  ,  .  »  -  B.  .  e.  .  ,  =  C 
i,l  i-l,l  i,l  i,l-l 


(3.35) 


The  condition  expressed  by  Equation  (3.35)  applies  when  i  and  j 
become  large  and  grid  points  are  far  from  boundaries.  Hence,  using 
Equations  (3.34)  and  (3.35),  with  the  assumption  that  the  influence  of 
the  boundary  points  is  neglected,  and  by  letting 

=  £,m  A^  ^  exp[i  (yiriAx  +  3tt  j  Az)  ]  (3.36) 

where  i2  =  -1  ,  we  can  obtain  the  decay  factor  £  per  iteration  for 
the  error  component  corresponding  to  y  and  8  as 


[0.5  C(l-a) K  ♦  CaA]  ♦  2CR 
*  ~  [0.5  C(l-a) K  +  CaA  +  S]  +  2CR 

where 

K  *  cos(ynAx)  cos(BtiAz)  , 


(3.37) 


„  .  5  YTtAX  .  ?  /0  7T  AZ 

A  =  2  sin2  sin^ 
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R  =  sin 


Y7TAX 

2 


sin 


BttAz 


cos 


ynAx 

2 


cos 


$7T  A  z 
2 


S  =  -  D  sin2  [  I^2L 


-  B  sin2 


BirAz 


In  order  to  keep  |c|  <  1  ,  the  iteration  parameter  a  has  to  be 
defined.  For  solving  self-adjoint  elliptic  difference  equations  Dupont, 
Kendall,  and  Rachford  (1968)  used  an  approximate  factorization  procedure. 
In  their  analysis,  a  Chebyshev  sequence  of  parameters  is  used  to  produce 
a  rapid  rate  of  convergence.  However,  neither  Equation  (2.15)  nor 
Equation  (2.16)  is  self-adjoint  and  their  method  for  determining  the 
iteration  parameter  a  is  not  applicable.  As  noted  by  Stone,  a  has 
to  be  less  than  1  and  the  repeated  use  of  any  a  in  the  vicinity  of  1 
results  in  |c|  >_  In  this  numerical  experiment,  both  a  transient 

heat  conduction  difference  equation  and  the  difference  equations  ex¬ 
pressed  in  Equations  (2.15),  (2.16),  and  (2.17)  have  been  investigated 
and  we  arrived  at  the  same  conclusions.  Stone  suggested  that  the 
individual  parameters  should  be  geometrically  spaced  in  the  manner 


1-a . 


=  (1-a  ) 

max7 


1/ CL- 1) 


l  =  0, _ L-l  , 


where  L  is  the  number  of  parameters  in  a  cycle  and  a  is  the 

luclX 

maximum  iteration  parameter.  Me  also  suggested  a  way  to  compute  amax: 


1 


a 

max 


rain 


2Ax2 


L  1 


KZAX2 

KXAz2 


1 


2Az2 

+  KXAz?~ 
KZAX2 


(3.38) 


where  KX  and  KZ  are  the  thermal  conductivities  in  the  x  and  z 
directions  in  the  thermal  conduction  equation 
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_a_ 

8x 


I 


+ 


9z 


Q 


(3.39) 


in  which  Q  is  the  strength  of  local  heat  source.  However,  neither 
Equation  (2.11)  nor  Equation  (2.12)  can  be  represented  in  as  simple  a 
form  as  Equation  (3.39),  and  even  if  it  could  be  done,  amax  calculated 
according  to  Equation  (3.38)  may  have  a  value  greater  than  1.  Hence 
in  this  study  we  simply  use  the  formula 


i  -  a  =  min 
max 


2Ax2 


2Az^ 


1  + 


Ax 

Az" 


1  + 


Az 

Ax^ 


(3.40) 


Experience  showed  that  a  ,  selected  by  using  Equation  (3.40),  has  no 
undesirable  effects.  Stone  also  indicated  a  possible  way  to  increase 
the  rate  of  convergence  by  use  of  double  alternate  iterations.  This 
process  can  be  carried  out  by  iterating  from  the  left- lower  corner  of 
the  field  to  its  right-upper  corner  and  then  from  the  left-upper  corner 
to  the  right-lower  corner.  However,  our  experience  indicated  that  the 
double  iteration  process  did  not  increase  the  rate  of  convergence  in 
the  problem  under  study,  and,  on  the  contrary,  it  may  sometimes  result 
in  a  decrease  of  the  convergence  rate.  For  this  reason,  we  apply  the 
process,  iterating  from  the  left-lower  corner  of  the  field  to  its  right- 
upper  corner.  The  rate  of  convergence  indicated  by  Equation  (3.37) 
applies  to  this  process. 

Before  the  iteration  method  described  above  was  applied  to  solve 
Equations  (2.15),  (2.16)  and  (2.17),  it  was  tried  on  a  test  problem  of 
a  flow  of  homogeneous,  viscous  fluid  with  governing  equations 
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+  3ji  iL  3jp  D£  3^£  3 2c 

3t  3z  3x  ”  3x  3z  "  .  2  +  .  2  * 

3x  3z 


and 


3  2  ip  3  2  tp 

„  2  +  „  2 

3x  3z 


C  , 


and  with  the  boundary  conditions 


(3.41) 


(3.42) 


rp  =  0  ,  iji  »  0  on  x  =  0.5,0<_z<0.5;  and  0  £  x  £  0.5  ,  z  =  0.5  , 
3;/3x  =  0  ,  3ip/3x  =  0  on  x  =  0  ,  0  <  z  <_  0.5  ,  (3.43) 

and 

3i;/3z  ,  3ij;/3z  =  0  on  0  <  x  <  0.5  ,  z  =  0 

An  exact  solution  of  Eqs.  (3.41)  and  (3.42)  with  boundary  conditions 
(3. 43) is  given  by 

ip  =  exp  (— 2lT2t )  COS7TX  COSTTZ  .  (3.44) 

The  nonlinear  terms  in  Equation  (3.41)  vanish  identically  for  the  exact 
solution  but,  in  the  computation  by  means  of  finite  difference  approxi¬ 
mation,  the  nonlinear  terms  do  not  quite  vanish;  thus,  any  instability 
caused  by  their  presence  should  be  detected. 

Since  Equation  (3.41)  is  identical  to  Equation  (2.1)  when  =  0 
and  Rg  =  1  ,  the  forward-backward  scheme  proposed  for  high  Reynolds 
number  was  net  used  because  of  the  low  Reynolds  number.  Instead,  the 
central  difference  form  is  used  for  3c/3x  and  3^/3z  .  For  the  numeri¬ 
cal  calculation,  Equation  (3.41)  is  expressed  in  an  implicit  finite 
difference  form  by 


35 


(3.45) 

where  the  superscripts  n  and  n+1  denote  the  time  steps  nAt  and 
(n«-l)At  ,  and  Equation  (3.42)  is  expressed  in  central  difference  form 
by 


(3.46) 

Figure  4  shows  grid  points  and  meshes  for  tho  test  problem.  In  the 
course  of  calculation,  we  used 

\p  =  COSTTX  costtz 
and 

c  =  TT2  COSIT X  COSTtZ 


as  the  initial  conditions  which  are  in  fact  the  exact  solution  presentt  ' 
in  Equation  (3.44)  when  t  =  0  .  The  calculation  field  contains  20  x  20 
grid  points  (including  the  boundary  points);  equal  mesh  size  is  used  for 
both  Ax  and  Az  ,  i.e..  Ax  =  Az  =  0.5/19;  the  time  increment  At  is 
0.002,  which  is  about  twelve  times  the  limit. 


(3.47) 
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given  by  stability  analysis  if  Equation  (3.41)  is  expressed  in  an 
explicit  finite  difference  form  and  if  it  is  assumed  that  the  stability 
criterion  is  determined  by  the  diffusion  terms  rather  than  by  the  con¬ 
vective  terms  of  Equation  (3.41).  In  fact,  it  can  be  shown  that  the 
required  time  increment  to  satisfy  the  stability  criterion  due  to  the 
presence  of  convective  terms  is 


Atd 


W/3zl 

Ax 


1  _ 

V 1 3i/.3.x- 

Az 


(3.48) 


which  is  much  greater  than  At^  expressed  in  Eq  tion  (3.47)  for  our 
present  problem. 

As  noted  in  section  3.3,  the  solution  expressed  in  Equation  (3.19) 
together  with  Equations  (3.20),  (3.21),  and  (3.22)  is  an  approximate 
solution  to  Equation  (3.T)  because  of  the  approximate  expressions 
Equation  (3.15)  and  Equation  (3.16)  for  T\  1  ^  and  T\+1  ^  j. 

Hence,  in  order  to  detect  the  effects  of  this  approximation  on  the  cal¬ 
culation,  we  calculate  Equation  (3.45)  once  for  a  given  time  step  and 
iteratively  calculate  Equation  (3.46)  twice.  From  this  process,  one 
can  also  detect  the  effects  of  the  presence  of  nonlinear  terms  on  the 
calculation.  Thus,  one  can  deduce  under  what  conditions  Equation  (3.45) 
has  to  be  calculated  iteratively  in  order  that  the  effects  due  to  the 
nonlinear  terms  may  be  taken  into  account.  As  a  summary  of  the  detailed 
calculation,  the  calculated  vorticity  and  stream  function  and  the  exact 
ones  at  the  point  x  =  1/19  and  z  =  1/19  are  listed  in  the  following 
table.  The  computer  used  was  a  CDC  6400. 


t 

0 

^exact 

9.60222368 

^calc. 

^exact 

0.97290862 

^calc.  1 

vcalc.  2 

0.002 

9.23052787 

9.22078202 

0.93524801 

0.92459860 

0.93348617 

0.004 

8.87322061 

8.85059308 

0.89904521 

0.88207465 

0.89408518 

0.006 

8.52974445 

8.49137856 

0.86424380 

0.84222221 

0.85587786 

0.008 

8.19956402 

8.14333278 

0.83078953 

0.80469707 

0.81901579 

0.01 

7.88216465 

7.80672216 

0.79863025 

0.76921420 

0.78360202 

0.012 

7.57705158 

7.48175327 

0.76771583 

0.73555790 

0.74967254 

0.014 

7.28374922 

7.16853732 

0.73799809 

0.70357176 

0.71722006 

0.016 

7.00180039 

6.86708870 

0.70943070 

0.67313486 

0.68621320 

0.018 

6.73076560 

6.57733479 

0.68196914 

0.64414754 

0.65660749 

0.02 

6.47022238 

6.29912885 

0.65557059 

0.61652354 

0.62835160 

0.022 

6.21976461 

6.03226294 

0.63019391 

0.59018586 

0.60139091 

In  the  table  4 1  1  .  and  iii  .  -  denote  the  stream  function  calcu- 

caic.  1  caic. 

lated  at  the  first  and  second  iteration,  respectively. 

Figure  5  presents  the  errors  in  calculated  stream  function  and 
vorticity  compared  to  their  exact  values  at  the  corresponding  times. 

At  t  =  0.002,  first  the  vorticity  was  calculated  by  Equation  (3.45) 
and  the’  this  newly  calculated  vorticity  was  used  in  Equation  (3.46) 
to  calculate  the  stream  function  iteratively  as  noted  by  4>  .  . 

Cel  X  0  •  X 

and  |l'calc  2*  These  calculations  indicate  that  the  result  of  even 
the  first  iteration  is  quite  close  to  the  corresponding  exact  solution. 
However,  because  of  the  presence  of  the  nonlinear  terms  the  error  grows 
monotonically  with  time  and,  thus,  it  is  concluded  that  for  solution  of 
nonlinear  difference  equations  such  as  Equation  (2.15),  or  Equation  (2.16) 
it  is  necessary  to  use  an  iterative  process  even  when  the  Reynolds  number 


is  small 
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Chipter  IV 

COMPUTATIONAL  PROCEDURE 


In  Chapter  III,  an  iterative  method  was  introduced  to  solve 
Equations  (2.15),  (2.16)  and  (2.17)  expressed  in  an  implicit  form.  For 
a  better  presentation  these  equations  are  rewritten  in  the  general  form 
expressed  in  Equation  (3.2): 


(i)  Temperature  field: 
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(ii)  Vorticity  field: 
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(iii)  Stream  function  field: 


L  d) ,  ♦  "L  d) .  -  2  ■  J~  ♦  d>  ♦  ■  1  ill 

>*’  l.J-1  ^  1-l.j  &x2  i*l.j 


♦  it ,  .  . 

a.?  *.’♦» 


1 


(4.3) 


40 


In  the  above  expressions,  the  superscripts  inside  parentheses  denote 
the  iteration  step  within  a  time  step  for  temperature  and  vorticity 
fields. 

In  the  course  of  computation  at  time  equal  to  zero,  a  uniform  flow 
starts  impulsively.  The  stream  function  field  is  calculation  by  letting 
C.  •  =  0  in  Equation  (4.3).  Of  course,  the  calculation  proceeds  itera- 

1  I  J 

tively  according  to  the  iterative  method  described  in  Chapter  III.  The 
initial  temperature  field  is  assumed  to  have  the  same  distribution  as 
the  initial  stream  function  field.  A  flow  chert  for  the  initial  solu¬ 
tion  is  presented  in  Figure  6.  From  the  initial  stream  function  initial 

velocities  u.  .  and  w.  .  can  be  calculated  to  determine  the  values 
i.J  i.J 

of  ,  Aj  ,  A^  and  A,.  according  to  the  expressions  Equation  (2.18). 
Subsequently,  the  temperature  field  at  the  interior  points  at  time  At 
can  be  computed  by  Equatio-  (4.1)  provided  the  boundary  conditions  in 
both  inflow  and  outflow  regions,  and  the  velocities  of  the  interior 
points  are  kept  at  the  initial  values;  then,  the  temperatures  on  the 
inflow  at”1  outflow  boundaries  are  extrapolated  by  expressions  Equations 
(2.30c)  and  (,2.31c).  The  vorticivy  field  at  the  interior  points  a.id  on 
inflow  and  outflow  boundaries  at  time  At  can  be  calculated  in  a  sim.- 
lar  manner.  However,  the  values  of  vorticity  on  the  obstacle  must  be 
evaluated  at  tine  level  At  by  Equations  (2.21),  (2.22)  and  (2.26) 
before  those  at  thr  interior  points  at  time  level  At  are  calculated. 
Subsequently  with  the  calculated  vorticity  at  time  At  ,  the  stream 
function  at  interior  points  at  time  At  is  iteratively  calculated  Hv 
Equation  (.4.3) .  Note  that  the  stream  function  n  the  inflow  and  outflow 
boundaries  is  that  from  the  (1-1)  iteration  when  the  ith  ite^:.c*on 
for  interior  points  is  processed  and  then  the  values  on  the  inflow  and 
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outflow  boundaries  at  the  £th  step  are  extrapolated  by  Equations  (2.30a) 
and  (2.31a)  from  those  at  interior  points  calculated  at  the  £th  iteration. 
The  iterative  process  continues  until  the  convergence  criterion 

Maxl*!1^  -  I  <  f.  ,  or  Max  |  1,3  <  6.  is  satisfied. 

1  1,3  1,3  1  (£+1)  1 

*i,j 

Mien  the  Reynolds  number  is  small  or,  equivalently,  when  the  nonlinear 

effects  are  not  important,  the  temperature,  vorticity,  and  stream  function 

at  time  2At  may  be  calculated  by  the  process  described  above.  However, 

usually  the  boundary  conditions  and  the  nonlinear  effects  are  so  closely 

coupled  that  an  iterative  process  must  be  used  for  the  temperature  and 

vorticity  fields.  To  carry  out  the  iterative  process,  we  register  the 

temperature  and  vorticity  fields,  calculated  as  described  above,  as 

T1’^  and  and  denote  the  associated  velocity  field  by  u?*^, 

i»3  1» J 


and  .  Subsequently,  T f  *  fm+  * ^  and  and  the  associated 
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velocity  field  and  w?’^  will  be  calculated  until  the  con- 
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are  satisfied.  At  this  stage  the  iterative  process  in  this  time  inter 
val  t  «  At  is  complete,  and  the  calculation  can  be  advanced  to  the 
next  time  interval  t  *  2At.  In  general,  the  iterative  process  can  he 
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applied  to  the  calculation  from  any  time  step,  say  the  nth  step,  to 
the  next  time  step,  say  the  n+lst  step.  Figure  6  shows  a  flow  chart 
describing  the  iterative  process  to  produce  the  initial  solution.  A 
general  flow  chart  for  a  transient  solution  from  the  nth  time  step  to 
the  n+lst  time  step  is  depicted  in  Figure  7. 
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Chapter  V 

RESULTS  AND  DISCUSSION 

Before  the  numerical  calculation  starts,  the  time  increment  and 
mesh  sizes  must  be  determined  from  different  points  of  view.  In  theory, 
for  an  implicit  scheme  there  is  no  restriction  on  the  time  increment. 
However,  because  of  boundary  conditions  and  of  the  nonlinearity  of  the 
equations,  an  iterative  process  is  always  needed  for  the  calculation 
of  vorticity  and  temperature.  If  a  large  time  increment  is  chosen,  it 
is  to  be  expected  that,  for  a  given  time  interval,  a  large  number  of 
iterations  will  be  required.  On  the  other  hand,  if  a  small  time  incre¬ 
ment  is  chosen  then,  for  a  given  time  interval,  fewer  iterations  will 
be  needed  but,  to  reach  an  asymptotic  state,  a  great  number  of  time  steps 
will  be  required.  There  is  no  way  to  determine  analytically  the  optimum 
time  increment.  The  only  accessible  way  is  to  conduct  a  numerical  experi¬ 
ment  and  to  find  out  the  optimum  time  increment  experimentally.  Similarly, 
the  optimum  mesh  size  can  also  be  determined  experimentally.  Of  course, 
the  computer  core  capacity  is  the  most  important  restriction  when  a  prob¬ 
lem  with  a  large  number  of  grid  points  is  calculated. 

Based  on  experience,  we  used  a  small  time  increment,  At  =  0.03  , 
for  the  first  calculation  starting  from  t  =  0,  since  the  flow  pattern 
changes  rapidly  right  after  the  fluid  starts  flowing.  After  that  the 
time  increment  At  was  increased  from  one  time  step  to  the  next  by  the 
formula 


At  =  At  ♦  0.02 


until  the  condition  At  >  0.2  is  reached. 
C.22,  was  used.  As  to  the  space  mesh  size, 


Thcreaf'or,  At  ■  0.2  or 
we  used  Ax  *  0.25  and 
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Az  =  0.75.  Smaller  mesh  sizes  have  been  tried,  but  the  limitation  in 
computer  storage  requires  the  computation  field  to  be  smaller  in  space 
and,  then,  a  lot  of  information  beyond  the  computation  field  is  lost. 

On  the  other  hand,  we  tried  to  use  larger  mesh  size,  but  the  increase 
in  truncation  errors  was  too  great  and  overshadowed  the  advantage  of  a 
larger  computation  field. 

Figures  8  to  14  are  the  computational  results  of  stream  function, 
vorticity,  and  temperature  at  different  time  steps  for  a  stratified  flow 
with  Reynolds  number  397,  Prandtl  number  10,  and  Richardson  number  1.58. 
The  figures  are  direct  prints  from  microfilm  plots  generated  by  CDC 
6600  computer.  These  dimensionless  numbers  are  defined  in  the  expres¬ 
sions  (2.10).  The  horizontal  and  vertical  coordinates  are  normalized 
by  the  height,  L  ,  of  the  vertical  fence  from  the  line  56  (Fig.  3); 
the  uniform  velocity  of  the  incoming  flow  is  the  characteristic  velocity; 
the  characteristic  temperature  difference  AT  is  defined  by  AT  =  dT/dz  L 
where  dT/dz  is  the  temperature  gradient  far  upstream  and  is  a  constant 
in  the  present  study.  Figure  8  shows  the  normalized  horizontal  and 
vertical  coordinates,  the  same  scale  being  used  for  both  coordinates. 

As  noted  previously,  the  fluid  starts  impulsively  as  an  irrotational 
flow,  and  thus,  at  the  early  stage  of  the  flow  development,  the  flow 
pattern  deviates  little  from  the  irrotational  flow  pattern.  Figure  8 
showing  the  stream  function  at  t  *  0,03  indicates  this  situation  very 
clearly  from  the  fact  that  the  streamlines  are  almost  symmetric  about 
the  vertical  axis  of  the  fence.  In  Figure  9  the  flow  pattern  at  t  ■ 

1.82  is  shown.  On  the  upstream  side  of  the  fence,  the  streamlines  are 
lifted  up,  while,  in  the  lee  side,  the  streamlines  converge  to  produce 
a  strong  downslope  current;  the  constant  temperature  lines  show  a 
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similar  pattern.  The  most  interesting  thing  is  the  development  of  vor- 
ticity.  In  the  case  of  homogeneous  flow  of  viscous  fluid,  vorticity  is 
confined  to  the  region  of  the  boundary  layer  and  it  may  diffuse  toward 
upstream  only  when  the  Reynolds  number  is  very  small.  However,  from 
Figure  9,  we  can  see  that  vorticity  is  not  only  diffused  in  the  wake 
region  but  also  exists  on  the  upstream  side  of  fence  as  a  result  of 
vorticity  creation  by  density  (or  temperature)  inhomogeneity. 

In  Figure  10,  lee-waves  are  being  generated  and  the  first  wave 
troughs  tilt  towards  upstream.  Figure  11  shows  a  clearer  flow  pattern 
of  lee-waves;  in  particular,  the  contour  lines  of  vorticity  present 
alternating  positive  and  negative  regions  indicating  the  troughs  and 
crests  of  lee-waves.  From  Figures  12  to  14,  one  can  see  the  development 
of  lee-waves  approaching  an  asymptotic  state.  In  Figure  14,  at  t  = 
15.68,  two  lee-waves  are  clearly  seen;  a  third  wave,  which  cannot  be 
observed  from  the  streamline  contours,  is  indicated  by  the  positive 
region  of  vorticity  in  the  left  hand  corner.  This  numerical  calcula¬ 
tion  can  be  compared  with  an  experimental  result  by  Pao  et  al.  (1968) 
for  a  cylinder  of  diameter  1.905  cm  roving  with  velocity  2.083  cm/sec 
in  stratified  salt  water  of  density  gradient  dp/dz  =  20.15  x  10-4  , 
and  it  will  be  found  that  the  flow  patterns  are  similar.  In  the  experi¬ 
ment  the  diameter  of  cylinder  is  the  characteristic  length,  while  in 
this  numerical  calculation,  the  fence  height  above  the  line  of  symmetry, 
corresponding  to  the  radius  of  cylinder  in  the  experiment,  serves  as 
the  characteristic  length.  For  further  comparison,  reference  may  be 
made  to  similar  flow  patterns  which  have  been  calculated  from  the  in- 
viscid  model  (Pao,  1969;  Miles,  1968). 
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Another  numerical  calculation  for  -  5000,  Pr  =  1,  and  = 

1.58  is  presented  in  Figures  15  to  21.  This  calculation  was  performed 
to  determine  whether  the  forward  and  backward  scheme  proposed  for  the 
vorticity  and  temperature  equations  of  stratified  flow  can  really  take 
care  of  high  Reynolds  numbers.  The  calculation  indicates  that  this 
scheme  gives  a  stable  computation  as  predicted.  In  the  series  cf 
patterns  in  Figures  15  to  19  with  the  corresponding  time  t  from  0.03 
to  18.42,  one  can  observe  a  similar  flow  development  to  that  described 
in  the  preceding  paragraphs  and,  of  course,  lee-waves  occur  downstream  of 
the  obstacle.  However,  in  Figures  20  and  21,  undesirable  disturbances 
appear  near  the  inflow  boundary.  In  order  to  understand  these  undesirable 
disturbances,  it  is  necessary  to  recall  the  numerical  experiments  on 
boundary  conditions. 

To  beg.’n  with,  we  used  a  set  of  prescribed  stream  function,  vorticity, 
and  temperature  on  the  inflow  boundary.  The  computational  field  contains 
180  x  30  grid  points  with  Ax  =  0.25,  and  At  =  0.25,  and  the  vertical 
fence  is  located  at  120  grid  points  from  the  inflow  boundary.  We  found 
that  this  kind  of  undesirable  disturbance  appeared  near  the  inflow 
boundary  as  early  as  t  about  4,  and  that  if  the  Richardson  number  was 
increased,  these  disturbances  might  appear  even  earlier.  At  first  it 
was  thought  that  these  disturbances  might  have  resulted  from  the  particu¬ 
lar  iterative  method.  In  order  to  clear  this  point,  a  flow  of  homogene¬ 
ous  fluid  over  the  same  obstacle  was  calculated  with  the  same  Reynolds 
number  and  boundary  conditions.  We  found  that,  for  this  homogeneous 
flow,  there  existed  no  such  undesirable  disturbance  near  the  upstream 
boundary  and  that  the  computation  could  be  carried  on  as  long  as  was 
desired,  i.c.,  the  computation  scheme  was  stable  and  the  iterative 
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method  was  proven  to  be  very  efficient  in  reducing  computational  efforts. 
Hence,  it  is  concluded  that  these  undesirable  disturbances  are  associated 
with  the  existence  of  stratification. 

Now,  one  may  ask  why  these  undesirable  disturbances  exist  near  the 
inflow  boundary  but  not  near  the  outflow  boundary.  Consider  a  transient 
flow  of  stratified  fluid  passing  over  an  obstacle.  The  disturbances 
caused  by  the  presence  of  this  obstacle  may  propagate  towards  both  up¬ 
stream  and  downstream  in  the  form  of  iternal  gravity  waves.  If  we  have 
suitable  boundary  conditions  for  both  inflow  and  outflow  boundaries, 
these  iternal  gravity  waves  in  the  computational  field  can  pass  through 
both  boundaries,  and  as  time  goes  on  the  asymptotic  state  can  be  achieved 
in  the  long  run.  However,  if  unsuitable  boundary  conditions  are  used, 
the  disturbances  propagating  from  the  obstacle  cannot  all  pass  through 
the  boundaries  but  may  reflect  back  from  the  boundaries  as  artificially 
introduced  disturbances,  lor  this  numerical  study,  the  predictor  formula 
expressed  in  Equation  (2.30)  and  Equation  (2.31)  was  used  to  extrapolate 
the  disturbances  from  the  interior  points  to  the  boundaries.  At  the 
outflow  boundary,  it  appears  that  the  outflow  can  efficiently  carry  the 
disturbance  away  from  the  computational  field.  However,  at  the  inflow 
boundary,  any  undesirable  disturbance  which  occurs  may  be  brought  back 
into  the  computational  field  by  the  inflow  as  an  artificial  disturbance 
which  will  finally  produce  undesirable  lee-waves  downstream  of  the 
inflow  boundary.  One  numerical  experiment  was  made  to  examine  the 
behavior  of  these  undesirable  lee-waves  and  it  was  found  that  they  could 
propagate  downstream  of  the  inflow  boundary.  Clearly,  this  phenomenon 
is  a  kind  of  computational  instability,  caused  by  the  interaction  of  the 
boundary  condition  and  the  flow  field  although,  based  on  the  stability 
analysis,  the  computational  scheme  is  stable  and  convergent. 
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In  short,  the  effects  of  boundary  conditions  on  the  computation  are 
dependent  on  the  flow  phenomena  themselves.  In  the  case  of  homogeneous 
fluid,  the  boundary  condition  at  the  inflow  boundary  usually  involves 
no  problem  since  disturbances  decay  exponentially  with  distance  but,  in 
the  case  of  stratified  fluid,  the  inflow  boundary  condition  becomes  cru¬ 
cial  since  disturbances  propagate  upstream  far  from  the  obstacle  and, 
therefore,  this  boundary  condition  must  receive  careful  treatment  if 
meaningful  results  are  to  be  obtained.  This  is  the  reason  why,  in  the 
problem  described  herein,  the  computational  region  upstream  of  the 
obstacle  contains  more  grid  points  than  that  downstream.  According  to 
experience  gained  during  this  numerical  experiment,  a  transformation  in 
coordinates  to  make  the  computational  field  extended  towards  far  upstream 
and  downstream  would  be  helpful.  Where  a  vertical  fence  as  an  obstacle 
is  concerned,  exponential  coordinates  and  elliptical  coordinates  may  be 
used  so  that  in  the  region  near  the  obstacle  a  fine  mesh  is  available 
and,  in  the  region  far  from  the  fence,  a  coarse  mesh  may  be  used.  Thus, 
not  only  better  accuracy  could  be  obtained  but  also  the  flow  boundaries 
of  the  computational  field  could  bo  pushed  as  far  away  as  the  computer 
storage  allows.  Because  of  limitation  in  computing  time,  these  special 
coordinates  have  not  been  tried.  However,  in  principle,  their  applica¬ 
tions  should  produce  no  problem  in  computation  if  the  iterative  method 
introduced  in  Chapter  III  is  used  since,  as  noted  previously,  this 
iterative  method  is  insensitive  to  the  coefficients  contained  in  Equa¬ 
tion  (3.1)  and,  in  fact,  the  transformation  in  coordinates  changes  only 
these  coefficients  but  not  the  basic  structure  of  Equation  (3.1). 

In  order  to  illustrate  how  the  computation  was  performed  for  both 
flow  conditions  described  above,  two  tables  are  listed  at  the  end  of 
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this  chapter.  Both  computations  were  made  for  square  meshes  in 
Cartesian  coordinates.  In  these  tables,  the  first  and  second  columns 
indicate  the  time  t  at  which  the  computation  was  made  and  the  time 
increment  At  ,  respectively.  The  third  column  is  the  tine  step  being 
computed  measured  from  the  instant  when  the  flow  starts  impulsively. 

The  fourth  and  fifth  columns  denote  the  iteration  numbers  for  tempera¬ 
ture  and  vorticity  and  for  stream  function  at  a  given  time  step,  while 
the  sixth  and  seventh  are  the  accumulated  iteration  numbers  for  the  cor¬ 
responding  variables.  The  first  computation  listed  in  Table  I  has  the 
following  convergence  criteria  for  iteration  of  temperature  and  vor¬ 
ticity  at  a  given  time  step  (n+1). 
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i,J  i.j  i.j 


and  for  stream  function  using  the  calculated  vorticity  at  the  (m+1) 
iteration  at  time  step  (n+1), 


.  I  ”  1  1 


Max 


i'lfJ  •  (Yvn- J-  i  i*|jl)i >  °-1  • 


i|>:  . 

i.j 


or 


Maxi^t1'*  -  <^0.001  for  *  0.1  . 

t.J  i.J  i»J“~ 
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For  the  second  computation  listed  in  Table  II,  the  convergence  criteria 
are: 


(a)  t  <  20 


Max 


Qn+1, (m+1)  Qn+1, (m) 

Jjj _ isj _ 

Qn+1, (m+1) 
i.j 


<  0.01  for 


n+1, (m+1) 


|0.  . 

1  i.J 


>  0.1 


I 


or 


Max |ev  . 
1  i.J 


n+ 1 , (m+ 1 ) 


0n+l,(m) 

i.j 


<  0.001 


for 


I  n+1, (m+1) 

1  i.j 


<  0.1 


» 


where  6.  .  =  T.  .  or  r.  .  , 

i.J  i.j  i.j 


and 


d’:1)  - 


Max 


"  J,'(Ul)1|,Ji  -°-01  f°r  "  0,1  * 


i.J 


or 


Maxi*}4!15  -  *}4J|  <  0.001  for  |*}4*1}|  <0.1  . 

A«J  1 .  J  1 .  J  — 


(bj  t  >  20 


»u+l* (m+1)  _  Qn+l,(m) 


j  I)1* - 1  i°-00‘  •  for  l»u‘t**1)l  *  »•> 


I.J 


or 


<  0.0001  for  Is”*!''"*1'!  <  0.1  . 

1  •  J  1  »  J  _  I.J  ~ 


where  tK  .  =  T.  .  or  r  .  , 

t.J  ».J  i.J 


and 
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Max 


* 


»JL 


* 


i' 


,(*♦!) 

i  J 


<  0.001  for  1 1^: 


(4+1) 


l.J 


>  0.’ 


or 


Max  i |i : 


(*♦1) 


i.J 


-  *> 


w 


i»r  - 


0.0001  for 


^+1> 

1,J 


<  0.1 


From  columns  4  and  5  of  Tables  I  and  II,  it  can  be  seen  that  when 
the  computation  starts,  a  number  of  iterations  is  needed  since,  in  the 
early  stages,  the  flow  is  highly  .time-dependent.  As  the  computation 
goes  through  the  middle  course,  the  iteration  numbers  are  as  low  as  2. 
Indeed  this  situation  clearly  indicates  that  the  iterative  method 
introduced  in  Chapter  III  is  really  a  good  one.  It  is  noted  that  in 
Pearsons  work  (1965),  a  smoothing  parameter  was  introduced  for  both  the 
boundary  and  interior  points  to  reduce  the  computational  efforts  by 
decreasing  iteration  numbers  for  vorticity  calculation  at  a  given  time 
step;  it  seems  to  us  that  the  iteration  parameter  a  for  the  present 
iteration  method  may  have  this  effect.  For  the  first  computation,  79 
time  steps  were  calculated  ar  l  30  minutes  of  CDC  6600  computer  time 
were  spent;  on  the  average,  each  time  step  took  25  seconds  approximately. 
One  thing  must  be  noted  about  the  last  row  of  the  Table  II  which  shows 
that  at  time  t  »  23. 4o  a  great  number  of  iterations  was  made  for  tempera¬ 
ture  and  vorticity  as  well  as  for  stream  function.  Comparison  of  these 
values  with  the  corresponding  ones  at  t  ■  20.52  indicates  that  the  great 
number  of  iterations  is  not  due  to  the  smaller  residuals  used  for  con¬ 
vergence  criteria  when  t  20.  In  fact,  one  may  relate  this  situation 
to  the  computational  instability  caused  by  the  unsuitable  boundary 
conditions  as  described  above.  In  other  words,  if  this  difficulty 
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particularly  involved  in  the  numerical  calculation  of  stratified  flows 
can  be  removed  by  some  special  coordinates,  the  iterative  method  is 
indeed  a  good  one-  for  the  numerical  calculation  of  general  problems  in 
geophysical  fluid  mechanics. 
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FIGURES 


Figure  I  Sketch  of  Grid  Points ,  and  Mesh  Sizes 


Flow 

Direction 


Figure  2  Sketch  of  Meoh  in  Vicinity  of  the 
Ve.  tical  Fence 
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Figure  3  Sketch  of  Interior  Field  and  Boundaries  of  Stratified  Flow  Problem 


t ,  time 


Figure  5  Error  Propagation  for  the  Test  Problem 


ENTRY 


Fig.  6.  Flow  Chart  for  Initial  Solution 
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Fig.  ; 


Flow  Chart  for  Tiansicnt  Solution 
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Patterns  At  t  =  20.52,  Re  =  397 
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